home *** CD-ROM | disk | FTP | other *** search
/ IRIX Base Documentation 2002 November / SGI IRIX Base Documentation 2002 November.iso / usr / share / catman / p_man / cat3 / SCSL / strsen.z / strsen
Encoding:
Text File  |  2002-10-03  |  11.6 KB  |  331 lines
  1.  
  2.  
  3.  
  4. SSSSTTTTRRRRSSSSEEEENNNN((((3333SSSS))))                                                          SSSSTTTTRRRRSSSSEEEENNNN((((3333SSSS))))
  5.  
  6.  
  7.  
  8. NNNNAAAAMMMMEEEE
  9.      STRSEN - reorder the real Schur factorization of a real matrix A =
  10.      Q*T*Q**T, so that a selected cluster of eigenvalues appears in the
  11.      leading diagonal blocks of the upper quasi-triangular matrix T,
  12.  
  13. SSSSYYYYNNNNOOOOPPPPSSSSIIIISSSS
  14.      SUBROUTINE STRSEN( JOB, COMPQ, SELECT, N, T, LDT, Q, LDQ, WR, WI, M, S,
  15.                         SEP, WORK, LWORK, IWORK, LIWORK, INFO )
  16.  
  17.          CHARACTER      COMPQ, JOB
  18.  
  19.          INTEGER        INFO, LDQ, LDT, LIWORK, LWORK, M, N
  20.  
  21.          REAL           S, SEP
  22.  
  23.          LOGICAL        SELECT( * )
  24.  
  25.          INTEGER        IWORK( * )
  26.  
  27.          REAL           Q( LDQ, * ), T( LDT, * ), WI( * ), WORK( * ), WR( * )
  28.  
  29. IIIIMMMMPPPPLLLLEEEEMMMMEEEENNNNTTTTAAAATTTTIIIIOOOONNNN
  30.      These routines are part of the SCSL Scientific Library and can be loaded
  31.      using either the -lscs or the -lscs_mp option.  The -lscs_mp option
  32.      directs the linker to use the multi-processor version of the library.
  33.  
  34.      When linking to SCSL with -lscs or -lscs_mp, the default integer size is
  35.      4 bytes (32 bits). Another version of SCSL is available in which integers
  36.      are 8 bytes (64 bits).  This version allows the user access to larger
  37.      memory sizes and helps when porting legacy Cray codes.  It can be loaded
  38.      by using the -lscs_i8 option or the -lscs_i8_mp option. A program may use
  39.      only one of the two versions; 4-byte integer and 8-byte integer library
  40.      calls cannot be mixed.
  41.  
  42. PPPPUUUURRRRPPPPOOOOSSSSEEEE
  43.      STRSEN reorders the real Schur factorization of a real matrix A =
  44.      Q*T*Q**T, so that a selected cluster of eigenvalues appears in the
  45.      leading diagonal blocks of the upper quasi-triangular matrix T, and the
  46.      leading columns of Q form an orthonormal basis of the corresponding right
  47.      invariant subspace.
  48.  
  49.      Optionally the routine computes the reciprocal condition numbers of the
  50.      cluster of eigenvalues and/or the invariant subspace.
  51.  
  52.      T must be in Schur canonical form (as returned by SHSEQR), that is, block
  53.      upper triangular with 1-by-1 and 2-by-2 diagonal blocks; each 2-by-2
  54.      diagonal block has its diagonal elemnts equal and its off-diagonal
  55.      elements of opposite sign.
  56.  
  57.  
  58.  
  59.  
  60.  
  61.  
  62.  
  63.                                                                         PPPPaaaaggggeeee 1111
  64.  
  65.  
  66.  
  67.  
  68.  
  69.  
  70. SSSSTTTTRRRRSSSSEEEENNNN((((3333SSSS))))                                                          SSSSTTTTRRRRSSSSEEEENNNN((((3333SSSS))))
  71.  
  72.  
  73.  
  74. AAAARRRRGGGGUUUUMMMMEEEENNNNTTTTSSSS
  75.      JOB     (input) CHARACTER*1
  76.              Specifies whether condition numbers are required for the cluster
  77.              of eigenvalues (S) or the invariant subspace (SEP):
  78.              = 'N': none;
  79.              = 'E': for eigenvalues only (S);
  80.              = 'V': for invariant subspace only (SEP);
  81.              = 'B': for both eigenvalues and invariant subspace (S and SEP).
  82.  
  83.      COMPQ   (input) CHARACTER*1
  84.              = 'V': update the matrix Q of Schur vectors;
  85.              = 'N': do not update Q.
  86.  
  87.      SELECT  (input) LOGICAL array, dimension (N)
  88.              SELECT specifies the eigenvalues in the selected cluster. To
  89.              select a real eigenvalue w(j), SELECT(j) must be set to w(j) and
  90.              w(j+1), corresponding to a 2-by-2 diagonal block, either
  91.              SELECT(j) or SELECT(j+1) or both must be set to either both
  92.              included in the cluster or both excluded.
  93.  
  94.      N       (input) INTEGER
  95.              The order of the matrix T. N >= 0.
  96.  
  97.      T       (input/output) REAL array, dimension (LDT,N)
  98.              On entry, the upper quasi-triangular matrix T, in Schur canonical
  99.              form.  On exit, T is overwritten by the reordered matrix T, again
  100.              in Schur canonical form, with the selected eigenvalues in the
  101.              leading diagonal blocks.
  102.  
  103.      LDT     (input) INTEGER
  104.              The leading dimension of the array T. LDT >= max(1,N).
  105.  
  106.      Q       (input/output) REAL array, dimension (LDQ,N)
  107.              On entry, if COMPQ = 'V', the matrix Q of Schur vectors.  On
  108.              exit, if COMPQ = 'V', Q has been postmultiplied by the orthogonal
  109.              transformation matrix which reorders T; the leading M columns of
  110.              Q form an orthonormal basis for the specified invariant subspace.
  111.              If COMPQ = 'N', Q is not referenced.
  112.  
  113.      LDQ     (input) INTEGER
  114.              The leading dimension of the array Q.  LDQ >= 1; and if COMPQ =
  115.              'V', LDQ >= N.
  116.  
  117.      WR      (output) REAL array, dimension (N)
  118.              WI      (output) REAL array, dimension (N) The real and imaginary
  119.              parts, respectively, of the reordered eigenvalues of T. The
  120.              eigenvalues are stored in the same order as on the diagonal of T,
  121.              with WR(i) = T(i,i) and, if T(i:i+1,i:i+1) is a 2-by-2 diagonal
  122.              block, WI(i) > 0 and WI(i+1) = -WI(i). Note that if a complex
  123.              eigenvalue is sufficiently ill-conditioned, then its value may
  124.              differ significantly from its value before reordering.
  125.  
  126.  
  127.  
  128.  
  129.                                                                         PPPPaaaaggggeeee 2222
  130.  
  131.  
  132.  
  133.  
  134.  
  135.  
  136. SSSSTTTTRRRRSSSSEEEENNNN((((3333SSSS))))                                                          SSSSTTTTRRRRSSSSEEEENNNN((((3333SSSS))))
  137.  
  138.  
  139.  
  140.      M       (output) INTEGER
  141.              The dimension of the specified invariant subspace.  0 < = M <= N.
  142.  
  143.      S       (output) REAL
  144.              If JOB = 'E' or 'B', S is a lower bound on the reciprocal
  145.              condition number for the selected cluster of eigenvalues.  S
  146.              cannot underestimate the true reciprocal condition number by more
  147.              than a factor of sqrt(N). If M = 0 or N, S = 1.  If JOB = 'N' or
  148.              'V', S is not referenced.
  149.  
  150.      SEP     (output) REAL
  151.              If JOB = 'V' or 'B', SEP is the estimated reciprocal condition
  152.              number of the specified invariant subspace. If M = 0 or N, SEP =
  153.              norm(T).  If JOB = 'N' or 'E', SEP is not referenced.
  154.  
  155.      WORK    (workspace/output) REAL array, dimension (LWORK)
  156.              On exit, if INFO = 0, WORK(1) returns the optimal LWORK.
  157.  
  158.      LWORK   (input) INTEGER
  159.              The dimension of the array WORK.  If JOB = 'N', LWORK >=
  160.              max(1,N); if JOB = 'E', LWORK >= M*(N-M); if JOB = 'V' or 'B',
  161.              LWORK >= 2*M*(N-M).
  162.  
  163.              If LWORK = -1, then a workspace query is assumed; the routine
  164.              only calculates the optimal size of the WORK array, returns this
  165.              value as the first entry of the WORK array, and no error message
  166.              related to LWORK is issued by XERBLA.
  167.  
  168.      IWORK   (workspace) INTEGER array, dimension (LIWORK)
  169.              IF JOB = 'N' or 'E', IWORK is not referenced.
  170.  
  171.      LIWORK  (input) INTEGER
  172.              The dimension of the array IWORK.  If JOB = 'N' or 'E', LIWORK >=
  173.              1; if JOB = 'V' or 'B', LIWORK >= M*(N-M).
  174.  
  175.              If LIWORK = -1, then a workspace query is assumed; the routine
  176.              only calculates the optimal size of the IWORK array, returns this
  177.              value as the first entry of the IWORK array, and no error message
  178.              related to LIWORK is issued by XERBLA.
  179.  
  180.      INFO    (output) INTEGER
  181.              = 0: successful exit
  182.              < 0: if INFO = -i, the i-th argument had an illegal value
  183.              = 1: reordering of T failed because some eigenvalues are too
  184.              close to separate (the problem is very ill-conditioned); T may
  185.              have been partially reordered, and WR and WI contain the
  186.              eigenvalues in the same order as in T; S and SEP (if requested)
  187.              are set to zero.
  188.  
  189.  
  190.  
  191.  
  192.  
  193.  
  194.  
  195.                                                                         PPPPaaaaggggeeee 3333
  196.  
  197.  
  198.  
  199.  
  200.  
  201.  
  202. SSSSTTTTRRRRSSSSEEEENNNN((((3333SSSS))))                                                          SSSSTTTTRRRRSSSSEEEENNNN((((3333SSSS))))
  203.  
  204.  
  205.  
  206. FURTHER DETAILS
  207.      STRSEN first collects the selected eigenvalues by computing an orthogonal
  208.      transformation Z to move them to the top left corner of T.  In other
  209.      words, the selected eigenvalues are the eigenvalues of T11 in:
  210.  
  211.                    Z'*T*Z = ( T11 T12 ) n1
  212.                             (  0  T22 ) n2
  213.                                n1  n2
  214.  
  215.      where N = n1+n2 and Z' means the transpose of Z. The first n1 columns of
  216.      Z span the specified invariant subspace of T.
  217.  
  218.      If T has been obtained from the real Schur factorization of a matrix A =
  219.      Q*T*Q', then the reordered real Schur factorization of A is given by A =
  220.      (Q*Z)*(Z'*T*Z)*(Q*Z)', and the first n1 columns of Q*Z span the
  221.      corresponding invariant subspace of A.
  222.  
  223.      The reciprocal condition number of the average of the eigenvalues of T11
  224.      may be returned in S. S lies between 0 (very badly conditioned) and 1
  225.      (very well conditioned). It is computed as follows. First we compute R so
  226.      that
  227.  
  228.                             P = ( I  R ) n1
  229.                                 ( 0  0 ) n2
  230.                                   n1 n2
  231.  
  232.      is the projector on the invariant subspace associated with T11.  R is the
  233.      solution of the Sylvester equation:
  234.  
  235.                            T11*R - R*T22 = T12.
  236.  
  237.      Let F-norm(M) denote the Frobenius-norm of M and 2-norm(M) denote the
  238.      two-norm of M. Then S is computed as the lower bound
  239.  
  240.                          (1 + F-norm(R)**2)**(-1/2)
  241.  
  242.      on the reciprocal of 2-norm(P), the true reciprocal condition number.  S
  243.      cannot underestimate 1 / 2-norm(P) by more than a factor of sqrt(N).
  244.  
  245.      An approximate error bound for the computed average of the eigenvalues of
  246.      T11 is
  247.  
  248.                             EPS * norm(T) / S
  249.  
  250.      where EPS is the machine precision.
  251.  
  252.      The reciprocal condition number of the right invariant subspace spanned
  253.      by the first n1 columns of Z (or of Q*Z) is returned in SEP.  SEP is
  254.      defined as the separation of T11 and T22:
  255.  
  256.                         sep( T11, T22 ) = sigma-min( C )
  257.  
  258.  
  259.  
  260.  
  261.                                                                         PPPPaaaaggggeeee 4444
  262.  
  263.  
  264.  
  265.  
  266.  
  267.  
  268. SSSSTTTTRRRRSSSSEEEENNNN((((3333SSSS))))                                                          SSSSTTTTRRRRSSSSEEEENNNN((((3333SSSS))))
  269.  
  270.  
  271.  
  272.      where sigma-min(C) is the smallest singular value of the
  273.      n1*n2-by-n1*n2 matrix
  274.  
  275.         C  = kprod( I(n2), T11 ) - kprod( transpose(T22), I(n1) )
  276.  
  277.      I(m) is an m by m identity matrix, and kprod denotes the Kronecker
  278.      product. We estimate sigma-min(C) by the reciprocal of an estimate of the
  279.      1-norm of inverse(C). The true reciprocal 1-norm of inverse(C) cannot
  280.      differ from sigma-min(C) by more than a factor of sqrt(n1*n2).
  281.  
  282.      When SEP is small, small changes in T can cause large changes in the
  283.      invariant subspace. An approximate bound on the maximum angular error in
  284.      the computed right invariant subspace is
  285.  
  286.                          EPS * norm(T) / SEP
  287.  
  288.  
  289. SSSSEEEEEEEE AAAALLLLSSSSOOOO
  290.      INTRO_LAPACK(3S), INTRO_SCSL(3S)
  291.  
  292.      This man page is available only online.
  293.  
  294.  
  295.  
  296.  
  297.  
  298.  
  299.  
  300.  
  301.  
  302.  
  303.  
  304.  
  305.  
  306.  
  307.  
  308.  
  309.  
  310.  
  311.  
  312.  
  313.  
  314.  
  315.  
  316.  
  317.  
  318.  
  319.  
  320.  
  321.  
  322.  
  323.  
  324.  
  325.  
  326.  
  327.                                                                         PPPPaaaaggggeeee 5555
  328.  
  329.  
  330.  
  331.